home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / shint.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  3.9 KB  |  136 lines
  1. /* specfunc/shint.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_expint.h>
  26.  
  27. #include "error.h"
  28.  
  29. #include "chebyshev.h"
  30. #include "cheb_eval.c"
  31.  
  32. /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  33.  
  34. /* based on SLATEC shi.f, W. Fullerton
  35.  
  36.  series for shi  on the interval  0.00000e+00 to  1.40625e-01
  37.                     with weighted error   4.67e-20
  38.                      log weighted error  19.33
  39.                    significant figures required  17.07
  40.                     decimal places required  19.75
  41. */
  42. static double shi_data[7] = {
  43.    0.0078372685688900950695,
  44.    0.0039227664934234563973,
  45.    0.0000041346787887617267,
  46.    0.0000000024707480372883,
  47.    0.0000000000009379295591,
  48.    0.0000000000000002451817,
  49.    0.0000000000000000000467
  50. };
  51. static cheb_series shi_cs = {
  52.   shi_data,
  53.   6,
  54.   -1, 1,
  55.   6
  56. };
  57.  
  58.  
  59. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  60.  
  61. int gsl_sf_Shi_e(const double x, gsl_sf_result * result)
  62. {
  63.   const double xsml = GSL_SQRT_DBL_EPSILON;  /* sqrt (r1mach(3)) */
  64.   const double ax   = fabs(x);
  65.  
  66.   if(ax < xsml) {
  67.     result->val = x;
  68.     result->err = 0.0;
  69.     return GSL_SUCCESS;
  70.   }
  71.   else if(ax <= 0.375) {
  72.     gsl_sf_result result_c;
  73.     cheb_eval_e(&shi_cs, 128.0*x*x/9.0-1.0, &result_c);
  74.     result->val  = x * (1.0 + result_c.val);
  75.     result->err  = x * result_c.err;
  76.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  77.     return GSL_SUCCESS;
  78.   }
  79.   else {
  80.     gsl_sf_result result_Ei;
  81.     gsl_sf_result result_E1;
  82.     int status_Ei = gsl_sf_expint_Ei_e(x, &result_Ei);
  83.     int status_E1 = gsl_sf_expint_E1_e(x, &result_E1);
  84.     result->val  = 0.5*(result_Ei.val + result_E1.val);
  85.     result->err  = 0.5*(result_Ei.err + result_E1.err);
  86.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  87.     if(status_Ei == GSL_EUNDRFLW && status_E1 == GSL_EUNDRFLW) {
  88.       GSL_ERROR ("underflow", GSL_EUNDRFLW);
  89.     }
  90.     else if(status_Ei == GSL_EOVRFLW || status_E1 == GSL_EOVRFLW) {
  91.       GSL_ERROR ("overflow", GSL_EOVRFLW);
  92.     }
  93.     else {
  94.       return GSL_SUCCESS;
  95.     }
  96.   }
  97. }
  98.  
  99.  
  100. int gsl_sf_Chi_e(const double x, gsl_sf_result * result)
  101. {
  102.   gsl_sf_result result_Ei;
  103.   gsl_sf_result result_E1;
  104.   int status_Ei = gsl_sf_expint_Ei_e(x, &result_Ei);
  105.   int status_E1 = gsl_sf_expint_E1_e(x, &result_E1);
  106.   if(status_Ei == GSL_EDOM || status_E1 == GSL_EDOM) {
  107.     DOMAIN_ERROR(result);
  108.   }
  109.   else if(status_Ei == GSL_EUNDRFLW && status_E1 == GSL_EUNDRFLW) {
  110.     UNDERFLOW_ERROR(result);
  111.   }
  112.   else if(status_Ei == GSL_EOVRFLW || status_E1 == GSL_EOVRFLW) {
  113.     OVERFLOW_ERROR(result);
  114.   }
  115.   else {
  116.     result->val  = 0.5 * (result_Ei.val - result_E1.val);
  117.     result->err  = 0.5 * (result_Ei.err + result_E1.err);
  118.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  119.     return GSL_SUCCESS;
  120.   }
  121. }
  122.  
  123. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  124.  
  125. #include "eval.h"
  126.  
  127. double gsl_sf_Shi(const double x)
  128. {
  129.   EVAL_RESULT(gsl_sf_Shi_e(x, &result));
  130. }
  131.  
  132. double gsl_sf_Chi(const double x)
  133. {
  134.   EVAL_RESULT(gsl_sf_Chi_e(x, &result));
  135. }
  136.